{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Trial 4: graph coarsening\n",
    "\n",
    "* First Python implementation of the greedy Metis and Graclus coarsening algorithms.\n",
    "* Results comparison with a previously developed matlab implementation.\n",
    "* Results comparison with the newer version in the `coarsening` module."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "METIS COARSENING IMPLEMENTATION AS PROPOSED IN:\n",
    "An incremental reseeding strategy for clustering\n",
    "X Bresson, H Hu, T Laurent, A Szlam, J von Brecht\n",
    "arXiv preprint arXiv:1406.3837\n",
    "3 May 2016"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import scipy.io\n",
    "import scipy.sparse\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    # CHECK PYTHON RESULTS WITH MATLAB CODE\n",
    "    folder = os.path.join('..', 'data', 'metis_matlab.mat')\n",
    "    mat = scipy.io.loadmat(folder)\n",
    "    W = mat['W']\n",
    "    W = scipy.sparse.csr_matrix(W)\n",
    "    rid = mat['rid']-1\n",
    "    rid = rid.T\n",
    "    rid = rid.squeeze()\n",
    "    #print(type(W))\n",
    "    #print(type(rid))\n",
    "    print(W.shape)\n",
    "    print(W.nnz)\n",
    "    #print(rid.shape)\n",
    "\n",
    "else:\n",
    "    N = 533\n",
    "    #np.random.seed(0)\n",
    "    rid = np.random.permutation(range(N))\n",
    "    W = np.random.uniform(0.01, 0.99, size=(N,N))\n",
    "    mask = np.random.uniform(size=(N,N))\n",
    "    W[mask<0.99] = 0\n",
    "    W = scipy.sparse.csr_matrix(W)\n",
    "    print(W.nnz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# INPUT\n",
    "# W = symmetric sparse weight matrix\n",
    "# maxsize = the number of nodes for the coarsest graph\n",
    "# OUTPUT\n",
    "# graph{1}: original graph of size N_1\n",
    "# graph{2}: coarser graph of size N_2 < N_1\n",
    "# etc...\n",
    "# graph{k}: corsest graph of Size N_k <...< N_2 < N_1\n",
    "# parents{i} is a vector of size N_i with entries ranging from 1 to N_{i+1}\n",
    "# which indicate the parents in the coarser graph{i+1} \n",
    "# nd_sz{i} is a vector of size N_i that contains the size of the supernode in the graph{i}\n",
    "# NOTE\n",
    "# if \"graph\" is a cell of size k, then \"parents\" will be a cell of size k-1\n",
    "\n",
    "def metis_coarsening(W,maxsize,rid):\n",
    "    \n",
    "    N = W.shape[0]\n",
    "    print('Size of original graph=',N)\n",
    "    parents = []\n",
    "    degree = W.sum(axis=0) - W.diagonal()\n",
    "    graphs = []\n",
    "    graphs.append(W)\n",
    "    supernode_size = np.ones(N)\n",
    "    nd_sz = [supernode_size]\n",
    "    count = 0\n",
    "    \n",
    "    while N > maxsize:\n",
    "        \n",
    "        count = count + 1;\n",
    "        print('level=',count)\n",
    "        \n",
    "        # CHOOSE THE WEIGHTS FOR THE PAIRING\n",
    "        # weights = ones(N,1)       # metis weights\n",
    "        weights = degree            # graclus weights\n",
    "        # weights = supernode_size  # other possibility\n",
    "        weights = weights.T\n",
    "        weights = np.array(weights)\n",
    "        weights = weights.squeeze()\n",
    "        \n",
    "         # PAIR THE VERTICES AND CONSTRUCT THE ROOT VECTOR\n",
    "        idx_row,idx_col,val = scipy.sparse.find(W) \n",
    "        perm = np.argsort(idx_row)\n",
    "        rr = idx_row[perm]\n",
    "        cc = idx_col[perm]\n",
    "        vv = val[perm]\n",
    "        cluster_id = one_level_coarsening(rr,cc,vv,rid,weights) # rr is ordered  \n",
    "        parents.append(cluster_id)\n",
    "        \n",
    "        # TO DO\n",
    "        # COMPUTE THE SIZE OF THE SUPERNODES AND THEIR DEGREE \n",
    "        #supernode_size = full(   sparse(cluster_id,  ones(N,1) , supernode_size )     )\n",
    "        #print(cluster_id)\n",
    "        #print(supernode_size)\n",
    "        #nd_sz{count+1}=supernode_size;\n",
    "    \n",
    "         # COMPUTE THE EDGES WEIGHTS FOR THE NEW GRAPH\n",
    "        nrr = cluster_id[rr]\n",
    "        ncc = cluster_id[cc]\n",
    "        nvv = vv\n",
    "        Nnew = int(cluster_id.max()) + 1\n",
    "        print('Size of coarser graph=',Nnew)\n",
    "        W = scipy.sparse.csr_matrix((nvv,(nrr,ncc)),shape=(Nnew,Nnew))\n",
    "        # Add new graph to the list of all coarsened graphs\n",
    "        graphs.append(W)\n",
    "        N = W.shape[0]\n",
    "        \n",
    "        # COMPUTE THE DEGREE (OMIT OR NOT SELF LOOPS)\n",
    "        degree = W.sum(axis=0)\n",
    "        #degree = W.sum(axis=0) - W.diagonal()\n",
    "    \n",
    "        # CHOOSE THE ORDER IN WHICH VERTICES WILL BE VISTED AT THE NEXT PASS\n",
    "        #[~, rid]=sort(ss);     # arthur strategy\n",
    "        #[~, rid]=sort(supernode_size);    #  thomas strategy\n",
    "        #rid=randperm(N);                  #  metis/graclus strategy   \n",
    "        ss = W.sum(axis=0).T\n",
    "        rid = [i[0] for i in sorted(enumerate(ss), key=lambda x:x[1])] # [~, rid]=sort(ss);\n",
    "        \n",
    "        \n",
    "    # Remove all diagonal entries in similarity matrices\n",
    "    for i in range(len(graphs)): \n",
    "        csr_setdiag_val(graphs[i])\n",
    "        scipy.sparse.csr_matrix.eliminate_zeros(graphs[i])\n",
    "    \n",
    "        \n",
    "    return graphs,parents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#http://nbviewer.ipython.org/gist/Midnighter/9992103\n",
    "def csr_setdiag_val(csr, value=0):\n",
    "    \"\"\"Set all diagonal nonzero elements\n",
    "    (elements currently in the sparsity pattern)\n",
    "    to the given value. Useful to set to 0 mostly.\n",
    "    \"\"\"\n",
    "    if csr.format != \"csr\":\n",
    "        raise ValueError('Matrix given must be of CSR format.')\n",
    "    csr.sort_indices()\n",
    "    pointer = csr.indptr\n",
    "    indices = csr.indices\n",
    "    data = csr.data\n",
    "    for i in range(min(csr.shape)):\n",
    "        ind = indices[pointer[i]: pointer[i + 1]]\n",
    "        j =  ind.searchsorted(i)\n",
    "        # matrix has only elements up until diagonal (in row i)\n",
    "        if j == len(ind):\n",
    "            continue\n",
    "        j += pointer[i]\n",
    "        # in case matrix has only elements after diagonal (in row i)\n",
    "        if indices[j] == i:\n",
    "            data[j] = value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Coarsen a graph given by rr,cc,vv.  rr is assumed to be ordered\n",
    "def one_level_coarsening(rr,cc,vv,rid,weights):\n",
    "    \n",
    "    nnz = rr.shape[0]\n",
    "    N = rr[nnz-1]+1\n",
    "    #print(nnz,N)\n",
    "    \n",
    "    marked = np.zeros(N)\n",
    "    rowstart = np.zeros(N)\n",
    "    rowlength = np.zeros(N)\n",
    "    cluster_id = np.zeros(N)\n",
    "    \n",
    "    oldval = rr[0]\n",
    "    count = 0\n",
    "    clustercount = 0\n",
    "    \n",
    "    for ii in range(nnz):\n",
    "        rowlength[count] = rowlength[count] + 1\n",
    "        if rr[ii] > oldval:\n",
    "            oldval = rr[ii]\n",
    "            rowstart[count+1] = ii\n",
    "            count = count + 1\n",
    "            \n",
    "    for ii in range(N):\n",
    "        tid = rid[ii]\n",
    "        if marked[tid]==0.0:\n",
    "            wmax = 0.0\n",
    "            rs = rowstart[tid]\n",
    "            marked[tid] = 1.0\n",
    "            bestneighbor = -1\n",
    "            for jj in range(int(rowlength[tid])):\n",
    "                nid = cc[rs+jj]\n",
    "                tval = (1.0-marked[nid]) * vv[rs+jj] *  (1.0/weights[tid]+ 1.0/weights[nid])\n",
    "                if  tval > wmax:\n",
    "                    wmax = tval\n",
    "                    bestneighbor = nid\n",
    "                    \n",
    "            cluster_id[tid] = clustercount;\n",
    "            \n",
    "            if  bestneighbor > -1:\n",
    "                cluster_id[bestneighbor] = clustercount\n",
    "                marked[bestneighbor] = 1.0\n",
    "        \n",
    "            clustercount = clustercount + 1\n",
    "            \n",
    "    return cluster_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "maxsize = 200\n",
    "N = W.shape[0]\n",
    "#rid = np.random.permutation(range(N))\n",
    "#print(N)\n",
    "#print(rid[0:10])\n",
    "\n",
    "graphs,parents = metis_coarsening(W.copy(),maxsize,rid)\n",
    "#print(graph)\n",
    "#print(parents)\n",
    "\n",
    "\n",
    "# CHECK RESULTS WITH MATLAB CODE\n",
    "graph0 = graphs[0]\n",
    "print(graph0.shape)\n",
    "print(graph0[0,:])\n",
    "\n",
    "graph1 = graphs[1]\n",
    "print(graph1.shape)\n",
    "print(graph1[0,:])\n",
    "\n",
    "graph2 = graphs[2]\n",
    "print(graph2.shape)\n",
    "print(graph2[0,:])\n",
    "\n",
    "parents0 = parents[0]\n",
    "print(parents0.shape)\n",
    "print(parents0[0:10])\n",
    "\n",
    "parents1 = parents[1]\n",
    "print(parents1.shape)\n",
    "print(parents1[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "from lib import coarsening\n",
    "\n",
    "graphs, parents = coarsening.metis(W, 2, rid)\n",
    "\n",
    "for i,A in enumerate(graphs):\n",
    "    M, M = A.shape\n",
    "    A = A.tocoo()\n",
    "    A.setdiag(0)\n",
    "    A = A.tocsr()\n",
    "    A.eliminate_zeros()\n",
    "    graphs[i] = A\n",
    "    print('Layer {0}: M_{0} = {1} nodes, {2} edges'.format(i, M, A.nnz))\n",
    "\n",
    "# CHECK RESULTS WITH MATLAB CODE\n",
    "graph0 = graphs[0]\n",
    "print(graph0.shape)\n",
    "print(graph0[0,:])\n",
    "\n",
    "graph1 = graphs[1].tocsr()\n",
    "print(graph1.shape)\n",
    "print(graph1[0,:])\n",
    "\n",
    "graph2 = graphs[2].tocsr()\n",
    "print(graph2.shape)\n",
    "print(graph2[0,:])\n",
    "\n",
    "parents0 = parents[0]\n",
    "print(parents0.shape)\n",
    "print(parents0[0:10])\n",
    "\n",
    "parents1 = parents[1]\n",
    "print(parents1.shape)\n",
    "print(parents1[0:10])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Python results\n",
    "\n",
    "Size of original graph= 533\n",
    "level= 1\n",
    "Size of coarser graph= 279\n",
    "level= 2\n",
    "Size of coarser graph= 147\n",
    "(533, 533)\n",
    "  (0, 18)\t0.810464124165\n",
    "  (0, 59)\t0.349678536711\n",
    "  (0, 60)\t0.591336229831\n",
    "  (0, 83)\t0.388420442335\n",
    "  (0, 105)\t0.255134781894\n",
    "  (0, 210)\t0.656852096558\n",
    "  (0, 226)\t0.900257809833\n",
    "  (0, 299)\t0.065093756932\n",
    "  (0, 340)\t0.810464124165\n",
    "  (0, 407)\t0.431454676752\n",
    "(279, 279)\n",
    "  (0, 44)\t1.63660876872\n",
    "  (0, 58)\t2.42459126058\n",
    "  (0, 71)\t0.186153138092\n",
    "  (0, 115)\t1.99313658383\n",
    "  (0, 167)\t1.24818832639\n",
    "  (0, 168)\t2.95891026039\n",
    "  (0, 179)\t0.388420442335\n",
    "  (0, 240)\t0.431454676752\n",
    "(147, 147)\n",
    "  (0, 21)\t5.1886032791\n",
    "  (0, 85)\t1.08484314421\n",
    "  (0, 87)\t0.353738954483\n",
    "  (0, 127)\t0.186153138092\n",
    "  (0, 135)\t1.88273900708\n",
    "  (0, 141)\t0.255134781894\n",
    "(533,)\n",
    "[  57.  148.  184.  237.   93.   93.   47.   28.  133.   71.]\n",
    "(279,)\n",
    "[ 127.    4.   88.  128.   50.  120.   54.  123.  146.   26.]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Matlab results\n",
    "\n",
    "ans =\n",
    "\n",
    "   (1,19)      0.8105\n",
    "   (1,60)      0.3497\n",
    "   (1,61)      0.5913\n",
    "   (1,84)      0.3884\n",
    "   (1,106)     0.2551\n",
    "   (1,211)     0.6569\n",
    "   (1,227)     0.9003\n",
    "   (1,300)     0.0651\n",
    "   (1,341)     0.8105\n",
    "   (1,408)     0.4315\n",
    "\n",
    "\n",
    "ans =\n",
    "\n",
    "   (1,45)      1.6366\n",
    "   (1,59)      2.4246\n",
    "   (1,72)      0.1862\n",
    "   (1,116)     1.9931\n",
    "   (1,168)     1.2482\n",
    "   (1,169)     2.9589\n",
    "   (1,180)     0.3884\n",
    "   (1,241)     0.4315\n",
    "\n",
    "\n",
    "ans =\n",
    "\n",
    "   (1,22)      5.1886\n",
    "   (1,86)      1.0848\n",
    "   (1,88)      0.3537\n",
    "   (1,128)     0.1862\n",
    "   (1,136)     1.8827\n",
    "   (1,142)     0.2551\n",
    "\n",
    "\n",
    "ans =\n",
    "\n",
    "    58\n",
    "   149\n",
    "   185\n",
    "   238\n",
    "    94\n",
    "    94\n",
    "    48\n",
    "    29\n",
    "   134\n",
    "    72\n",
    "\n",
    "\n",
    "ans =\n",
    "\n",
    "   128\n",
    "     5\n",
    "    89\n",
    "   129\n",
    "    51\n",
    "   121\n",
    "    55\n",
    "   124\n",
    "   147\n",
    "    27"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
